# -*- coding: utf-8 -*-
"""
Created on Sun Apr 10 20:58:33 2022

@author: juju
"""

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt

data_csv=pd.read_csv("D:/03研究生学习/数学建模/示例数据/附件1-区域15分钟负荷数据.csv",usecols=[1])

# 数据预处理
data_csv = data_csv.dropna() # 滤除缺失数据
dataset = data_csv.values # 获得csv的值
dataset = dataset.astype('float32')
max_value = np.max(dataset) # 获得最大值
min_value = np.min(dataset) # 获得最小值
dataset = list(map(lambda x: (x-min_value) / (max_value-min_value), dataset)) # 标准化

plt.plot(dataset)
plt.show()

#设置X,Y数据集。以look_back=480(5天)为准，取第一个5天为数组，形成data_X,取第二个5天作为预测值，形成data_Y，完成训练集的提取。

def create_dataset(dataset, look_back=480):
    dataX, dataY = [], []
    for i in range(len(dataset) - look_back*2):
        a = dataset[i:(i + look_back)]
        dataX.append(a)
        dataY.append(dataset[i + look_back:i + look_back*2])
    return np.array(dataX), np.array(dataY)

# 创建好输入输出
data_X, data_Y = create_dataset(dataset)

# 划分训练集和测试集，70% 作为训练集
train_size = int(len(data_X) * 0.7)
test_size = len(data_X) - train_size
train_X = data_X[:train_size]
train_Y = data_Y[:train_size]
test_X = data_X[train_size:]
test_Y = data_Y[train_size:]

# 设置LSTM能识别的数据类型，形成tran_X的一维两个参数的数组，train_Y的一维一个参数的数组。并转化为tensor类型
import torch
# 把list转numpy三维数组，第一维自适应
train_X = train_X.reshape(-1, 1, 480)
train_Y = train_Y.reshape(-1, 1, 480)
test_X = test_X.reshape(-1, 1, 480)
# numpy数组转tensor
train_x = torch.from_numpy(train_X)
train_y = torch.from_numpy(train_Y)
test_x = torch.from_numpy(test_X)

# 建立LSTM模型，第一层为LSTM神经网络，第二层为一个全连接层。
from torch import nn
from torch.autograd import Variable

class lstm(nn.Module):
    def __init__(self,input_size=480,hidden_size=480,output_size=480,num_layer=2):
        super(lstm,self).__init__()
        self.layer1 = nn.LSTM(input_size,hidden_size,num_layer)
        self.layer2 = nn.Linear(hidden_size,output_size)

    def forward(self,x):
        x,_ = self.layer1(x)
        s,b,h = x.size()
        x = x.view(s*b,h)
        x = self.layer2(x)
        x = x.view(s,b,-1)
        return x
    
model = lstm(480, 480,480,2)
model = model.cuda()

# 设置交叉熵损失函数和自适应梯度下降算法
criterion = nn.MSELoss()
optimizer = torch.optim.Adam(model.parameters(), lr=1e-2)

# 开始训练
for e in range(2000):
    var_x = Variable(train_x).cuda()
    var_y = Variable(train_y).cuda()
    # 前向传播
    out = model(var_x)
    loss = criterion(out, var_y)
    # 反向传播
    optimizer.zero_grad()
    loss.backward()
    optimizer.step()

    if (e + 1) % 100 == 0: # 每 100 次输出结果
        print('Epoch: {}, Loss: {:.5f}'.format(e + 1, loss.item()))
        
model = model.eval() # 转换成测试模式
# data_X = data_X.reshape(-1, 1, 2)
# data_X = torch.from_numpy(data_X)
var_test_x = Variable(test_x).cuda()
pred_test_y = model(var_test_x) # 测试集的预测结果
# 改变输出的格式
# pred_test = pred_test.view(-1).data.numpy()
pred_test_Y = pred_test_y.view(-1).data.cpu().numpy().reshape(-1,480,1)

# 取最后一段5天预测的结果和实际对比，画出测试集中实际结果和预测的结果
pred_last_test_Y=pred_test_Y[-1]
last_test_Y=test_Y[-1]

plt.plot(pred_last_test_Y, 'r', label='prediction')
plt.plot(last_test_Y, 'b', label='real')
plt.legend(loc='best')
plt.show()

# 分析一下误差
# 均方误差
MSE = np.linalg.norm(last_test_Y-pred_last_test_Y, ord=2)**2/len(last_test_Y)
# 平均绝对误差
MAE = np.linalg.norm(last_test_Y-pred_last_test_Y, ord=1)/len(last_test_Y)
# 平均绝对百分比误差
MAPE = np.mean(np.abs((last_test_Y-pred_last_test_Y) / last_test_Y)) * 100
# 模型的准确率
Accuracy=(1-np.sqrt(np.mean(((last_test_Y-pred_last_test_Y) / last_test_Y)**2)))*100

print("MSE={}\nMAE={}\nMAPE={}%\nAccuracy={}%\n".format(MSE,MAE,MAPE,Accuracy))

# 使用滑动窗口预测未来的90天，也就是18个5天(480个15分钟)
pred_list=[]
x=dataset[-480:]
for i in range(18):
    input_x=[]
    input_x.append(x)
    input_x=np.array(input_x)
    input_x=input_x.reshape(-1,1,480)
    var_input_x=Variable(torch.from_numpy(input_x)).cuda()
    pred=model(var_input_x)
    pred=pred.view(-1).data.cpu().numpy()
    x=pred
    pred_list.append(pred)

prediction=np.array(pred_list).reshape(-1)

# 画出历史15天的图形和预测的90天的图形
plt.plot(range(len(dataset)),dataset, 'b', label='history')
plt.plot(range(len(dataset),len(dataset)+len(prediction)),prediction, 'r', label='prediction')
plt.legend(loc='best')
plt.show()

# 预测数据逆归一化
pred_power=prediction*(max_value-min_value)+min_value

def calc_date(index):
    if index>=0 and index<=16:
        return "2018/1/"+str(index+15)
    elif index<=44:
        return "2018/2/"+str(index-16)
    elif index<=75:
        return "2018/3/"+str(index-44)
    elif index<=105:
        return "2018/4/"+str(index-75)
    else:
        return "error"

# 一天96个数，一共90天
pred_power_array=np.split(pred_power,90)
pred_df= pd.DataFrame(columns=['date', 'max_power', 'max_time','min_power', 'min_time'])
for i in range(len(pred_power_array)):
    date=calc_date(i)
    max_power=np.max(pred_power_array[i])
    max_minute=np.argmax(pred_power_array[i])*15
    max_time="{}:{}".format(max_minute//60,max_minute%60)
    
    min_power=np.min(pred_power_array[i])
    min_minute=np.argmin(pred_power_array[i])*15
    min_time="{}:{}".format(min_minute//60,min_minute%60)
    pred_df.loc[len(pred_df.index)]=[date,max_power,max_time,min_power,min_time]
    
# 保存csv
pred_df.to_csv("./未来3个月预测.csv")
